//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2009 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : defines.cpp
//----------------------------------------------------------------------------


#include "matvec.h"


/// matrix vector product
Array2D<REAL> matvec(const std::vector<Array2D<REAL> > &A, const Array2D<REAL> &X)
{
    Array2D<REAL> Az(X);

    const int nx1 = X.dim1()-2;
    const int nx2 = X.dim2()-2;

    //std::cout << nx1 <<" " << nx2 <<std::endl;
    const Array2D<REAL> cc = A[0];
    const Array2D<REAL> cs = A[1];
    const Array2D<REAL> cw = A[2];
    const Array2D<REAL> cn = A[3];
    const Array2D<REAL> ce = A[4];


    /*for (int ii = 1; ii <= nx1; ++ii)
    {
        Az[ii][0] = cs[ii    ][1] * X[ii    ][1] +
                ce[ii ][0] * X[ii + 1][0] +
                cw[ii    ][0] * X[ii - 1][0] +
                cc[ii][0]*X[ii][0];
    }

    for (int jj = 1; jj <= nx2; ++jj)
    {
        Az[0][jj] = cn[0][jj] * X[0][jj + 1] +
                cw[1][jj    ] * X[1][jj    ] +
                cs[0][jj    ] * X[0][jj - 1] +
                X[0][jj]* cc[0][jj];
    }

    { //------------------------------------------------------------------------
        Az[0][0] = cc[0][0] * X [0][0] +
                cn[0][1] * X[0][1] +
                ce[1][0] * X[1][0];
    }*/
    int i1,j1;


#pragma omp for private(i1,j1)
    for (i1 = 1; i1 < nx1; ++i1)
    {
        int i1e = i1+1;
        int i1w = i1-1;
        for (j1 = 1; j1 < nx2; ++j1)
        {
            int j1n = j1+1;
            int j1s = j1-1;

            Az[i1][j1] = cn[i1    ][j1 ] * X[i1    ][j1n] +
                    ce[i1 ][j1    ] * X[i1e][j1    ] +
                    cs[i1    ][j1    ] * X[i1    ][j1s] +
                    cw[i1    ][j1    ] * X[i1w][j1    ]
                    + X[i1][j1]*cc[i1][j1];
        }
    }

    return Az;

}


// Adds the two vectors and returns the result
Array2D<REAL> MatAdd(const Array2D<REAL> &Vec1, const Array2D<REAL> &Vec2)
{
    Array2D<REAL> VecAdd(Vec1);
    const int nx1 = Vec1.dim1();
    const int nx2 = Vec1.dim2();


    if (Vec1.dim1()!=Vec2.dim1() && Vec1.dim2()!=Vec2.dim2() )
    {
        printf(" Matrix Addition- Dimensions do not agree");
        return VecAdd;
    }
    int i1,j1;



#pragma omp parallel for private(i1,j1)

    for (i1 = 0; i1 < nx1; i1++)
    {
        for (j1 = 0; j1 < nx2; j1++)
        {
            VecAdd[i1][j1] = Vec1[i1][j1] + Vec2[i1][j1];
        }
    }

    return VecAdd;
}


Array2D<REAL> MatSub(const Array2D<REAL> &Vec1, const Array2D<REAL> &Vec2)
{
    Array2D<REAL> VecSub(Vec1);
    const int nx1 = Vec1.dim1();
    const int nx2 = Vec2.dim2();

    if (Vec1.dim1()!=Vec2.dim1() && Vec1.dim2()!=Vec2.dim2() )
    {
        printf(" Dot Product- Dimensions do not agree");
        return VecSub;
    }
    int i1,j1;


#pragma omp parallel for private(i1,j1)

    for (i1 = 0; i1 < nx1; i1++)
    {
        for (j1 = 0; j1 < nx2; j1++)
        {

            VecSub[i1][j1] = Vec1[i1][j1] - Vec2[i1][j1];
        }
    }

    return VecSub;
}



Array2D<REAL> scalarmult(const REAL &mult, const Array2D<REAL> &Vec)
{
    Array2D<REAL> Vec_mult(Vec);

    int i1,j1;

#pragma omp parallel for private(i1,j1)

    for (i1 = 0; i1 < Vec.dim1(); i1++)
    {
        for (j1 = 0; j1 < Vec.dim2(); j1++)
        {

            Vec_mult[i1][j1] = Vec[i1][j1]*mult;
        }
    }

    // Vec_mult*=mult;
    return Vec_mult;
}

// Naive implementation of the Dot Product
REAL dotproduct(const Array2D<REAL> &Vec1, const Array2D<REAL> &Vec2)
{
    REAL sum =0;

    const int nx1 = Vec1.dim1()-2;
    const int nx2 = Vec1.dim2()-2;

    if (Vec1.dim1()!=Vec2.dim1() && Vec1.dim2()!=Vec2.dim2() )
    {
        printf(" Dot Product- Dimensions do not agree");
        return 0;
    }

    int i1,j1;

#pragma omp parallel for private(i1,j1) reduction (+:sum)
    for (i1 = 1; i1 < nx1; i1++)
    {
        for (j1 = 1; j1 < nx2; j1++)
        {

            sum += Vec1[i1][j1]*Vec2[i1][j1];
        }
    }
    return sum;
}


REAL norm(const Array2D<REAL> &Vec1)
{
    REAL Norm=0;
    Norm = sqrt(dotproduct(Vec1,Vec1));
    return Norm;
}

REAL norm(const Array2D<REAL> &Vec1,const Array2D<REAL> &Vec2)
{
    REAL Norm=0;
    Norm = dotproduct(Vec1,Vec1) + dotproduct(Vec2,Vec2);
    return sqrt(Norm);
}

void neumannBoundary(Array2D<REAL> * m_vec, int Nx1, int Nx2)
{
    Array2D<REAL> &m_Zeta = *m_vec;
    int nx1=m_Zeta.dim1()-2;
    int nx2=m_Zeta.dim2()-2;

    /* north and south*/

#pragma omp parallel for
    for (int i1 =1; i1< Nx1 ; ++ i1)
    {
        m_Zeta[i1][nx2] = m_Zeta[i1][Nx2];
        m_Zeta[i1][0] = m_Zeta[i1][1];

    }

    /* east and west*/
#pragma omp parallel for
    for (int i2 =1; i2< Nx2 ; ++ i2)
    {
        m_Zeta[nx1][i2] = m_Zeta[Nx1][i2];
        m_Zeta[0][i2] = m_Zeta[1][i2];
    }
}

void constructMatrix(int Nx1, int Nx2, REAL coeff_x, REAL coeff_y, REAL coeff_center, std::vector<Array2D<REAL> > &mat )
{
    Array2D<REAL> mat_c(Nx1+3, Nx2+3);
    Array2D<REAL> mat_w(Nx1+3, Nx2+3);
    Array2D<REAL> mat_s(Nx1+3, Nx2+3);
    mat_c=1.0;
    mat_w=0;
    mat_s=0;

    for (int i2 = 1; i2 <= Nx2; ++i2)
    {
        for (int i1 = 1; i1 <= Nx1; ++i1)
        {
            mat_c[i1][i2] = coeff_center;
            mat_w[i1][i2] = coeff_x;
            mat_s[i1][i2] = coeff_y;
        }
    }

    for (int i1 = 1; i1 <= Nx1; ++i1)
    {
        mat_c[i1][Nx2] +=mat_s[i1][Nx2];
        //mat_zz_cn[i1][Nx2]=0;
        mat_c[i1][1] +=mat_s[i1][1];
        mat_s[i1][1]=0;

    }

    /* east and west*/
    for (int i2 = 1; i2 <= Nx2; ++i2)
    {
        mat_c[Nx1][i2] +=mat_w[Nx1][i2];
        //  mat_zz_e[Nx1][i2]=0;
        mat_c[1][i2] +=mat_w[1][i2];
        mat_w[1][i2]=0;
    }


    Array2D<REAL> cn = mat_s;
    Array2D<REAL> ce = mat_w;

    for (int jj = 0; jj <= Nx2 + 1; ++jj)
    {
        for (int ii = 0; ii <= Nx1 + 1; ++ii)
        {
            cn[ii][jj] = mat_s[ii ][jj + 1]; // Make use of the fact that the matrix is Skew Symmetric.
            ce[ii][jj] = mat_w[ii + 1][jj ]; // Make use of the fact that the matrix is Skew Symmetric.
        }
    }

    mat.push_back(mat_c);
    mat.push_back(mat_s);
    mat.push_back(mat_w);
    mat.push_back(cn);
    mat.push_back(ce);

}



